data<-read.table("C:/Users/franc/Desktop/mammalia-macaque-dominance.edges", quote="\"", comment.char="")
mmd_graph <- graph_from_data_frame(data, directed=FALSE)
data1 <- data %>% transmute(macaque=V1, interactions=V2, weights=V3)
subjects <- data %>% transmute(macaque=V1, interactions=V2, weights=V3) %>%
bind_rows(data %>% transmute(macaque=V1, interactions=V2, weights=V3)) %>%
distinct() %>% arrange(macaque)
subjects
Here also, we’ve considered our data as an association network by connecting the nodes that belong to the ‘same’ group, so the individuals that assume the same behaviour may be considered to be ‘associates’ and therefore connected in a network. This is obviously an assumption in terms of dominance due to the lack of information related to macaques.
mmd_graph <- graph_from_data_frame(data1, directed=FALSE)
mmd_graph
## IGRAPH 2f4f722 UN-- 62 1167 --
## + attr: name (v/c), weights (e/n)
## + edges from 2f4f722 (vertex names):
## [1] 1--2 1--3 1--8 1--10 1--11 1--12 1--13 1--14 1--16 1--18 1--20 1--21
## [13] 1--22 1--23 1--24 1--26 1--30 1--31 1--33 1--35 1--36 1--37 1--38 1--40
## [25] 1--42 1--43 1--45 1--47 1--48 1--49 1--50 1--51 1--52 1--55 1--56 1--57
## [37] 1--58 1--59 2--10 2--12 2--24 2--29 2--30 2--35 2--36 2--37 2--38 2--40
## [49] 2--41 2--44 2--45 2--47 2--48 2--50 2--51 2--53 2--55 2--56 2--58 2--59
## [61] 2--61 3--5 3--6 3--8 3--10 3--11 3--13 3--14 3--16 3--17 3--19 3--20
## [73] 3--21 3--24 3--25 3--26 3--27 3--28 3--29 3--30 3--32 3--33 3--34 3--35
## [85] 3--36 3--37 3--38 3--39 3--41 3--42 3--44 3--45 3--46 3--48 3--49 3--51
## + ... omitted several edges
Firstly, let’s have a look at the structure of the graph:
str(mmd_graph)
## Class 'igraph' hidden list of 10
## $ : num 62
## $ : logi FALSE
## $ : num [1:1167] 1 2 7 9 10 11 12 13 15 17 ...
## $ : num [1:1167] 0 0 0 0 0 0 0 0 0 0 ...
## $ : num [1:1167] 0 1 61 104 62 133 2 63 105 134 ...
## $ : num [1:1167] 0 1 2 3 4 5 6 7 8 9 ...
## $ : num [1:63] 0 0 1 2 2 4 5 6 11 14 ...
## $ : num [1:63] 0 38 61 104 133 165 181 200 231 262 ...
## $ :List of 4
## ..$ : num [1:3] 1 0 1
## ..$ : Named list()
## ..$ :List of 1
## .. ..$ name: chr [1:62] "1" "2" "3" "4" ...
## ..$ :List of 1
## .. ..$ weights: int [1:1167] 3 1 1 2 1 2 1 5 1 1 ...
## $ :<environment: 0x000001daa52cf550>
In terms of response, the dataset contains information about dominance behavior among macaque monkeys, the response variable could be a column indicating each individual based on its dominant or submissive status.
response <- data1$macaque
response
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [25] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
## [49] 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3
## [73] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [97] 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [121] 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5
## [145] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6
## [169] 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7
## [193] 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
## [217] 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9
## [241] 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 10 10
## [265] 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
## [289] 10 10 10 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
## [313] 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 12 12 12 12 12 12 12
## [337] 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## [361] 12 12 12 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13 13 13 13 13 13 13
## [385] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
## [409] 13 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14
## [433] 14 14 14 14 14 14 14 14 14 14 14 14 14 14 15 15 15 15 15 15 15 15 15 15
## [457] 15 15 15 15 15 15 15 15 15 15 15 15 15 15 16 16 16 16 16 16 16 16 16 16
## [481] 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16
## [505] 16 16 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 18 18
## [529] 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 19 19
## [553] 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19
## [577] 19 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20
## [601] 20 20 20 20 20 20 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21
## [625] 21 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22 22
## [649] 22 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23
## [673] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
## [697] 24 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25
## [721] 25 25 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26 26
## [745] 26 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27
## [769] 27 28 28 28 28 28 28 28 28 28 28 28 28 28 28 28 28 28 28 29 29 29 29 29
## [793] 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 30 30 30 30 30 30 30 30 30
## [817] 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 31 31 31 31 31
## [841] 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 32 32 32 32 32 32
## [865] 32 32 32 32 32 32 32 32 32 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33
## [889] 33 33 34 34 34 34 34 34 34 34 34 34 34 34 34 34 34 34 34 34 35 35 35 35
## [913] 35 35 35 35 35 35 35 35 35 35 35 35 35 35 35 35 35 36 36 36 36 36 36 36
## [937] 36 36 36 36 36 36 36 36 36 36 36 37 37 37 37 37 37 37 37 37 37 37 37 37
## [961] 37 37 37 38 38 38 38 38 38 38 38 38 38 38 38 39 39 39 39 39 39 39 39 39
## [985] 39 39 39 39 39 39 39 39 39 39 40 40 40 40 40 40 40 40 40 40 40 40 40 40
## [1009] 40 40 40 41 41 41 41 41 41 41 41 41 41 41 41 41 41 41 41 41 41 42 42 42
## [1033] 42 42 42 42 42 42 42 42 42 42 42 43 43 43 43 43 43 43 43 43 43 44 44 44
## [1057] 44 44 44 44 44 44 44 44 45 45 45 45 45 45 45 45 45 45 45 46 46 46 46 46
## [1081] 46 46 46 46 46 47 47 47 47 47 47 47 47 47 47 48 48 48 48 48 48 48 48 48
## [1105] 48 48 48 49 49 49 49 49 49 49 49 49 49 50 50 50 51 51 51 51 51 52 52 52
## [1129] 52 52 53 53 53 53 53 53 53 54 54 54 54 54 55 55 55 55 55 55 55 56 56 56
## [1153] 56 56 56 57 57 57 57 58 58 58 59 59 59 60 60
The output obtained for the response variable is a sequence of numbers from 1 to 40 (representative for the whole 62 involved) where each number represents a individual, which can be considered as the response variable in the dataset. So in this case, the response variable is related to the different levels of dominance within the macaque population.
Then let’s explore the attributes of the mmd_graph object using the V(mmd_graph) function for nodes and the E(mmd_graph) for edges in order to investigate which attribute might serve as the response variable.
V(mmd_graph)
## + 62/62 vertices, named, from 2f4f722:
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62
E(mmd_graph)
## + 1167/1167 edges from 2f4f722 (vertex names):
## [1] 1--2 1--3 1--8 1--10 1--11 1--12 1--13 1--14 1--16 1--18 1--20 1--21
## [13] 1--22 1--23 1--24 1--26 1--30 1--31 1--33 1--35 1--36 1--37 1--38 1--40
## [25] 1--42 1--43 1--45 1--47 1--48 1--49 1--50 1--51 1--52 1--55 1--56 1--57
## [37] 1--58 1--59 2--10 2--12 2--24 2--29 2--30 2--35 2--36 2--37 2--38 2--40
## [49] 2--41 2--44 2--45 2--47 2--48 2--50 2--51 2--53 2--55 2--56 2--58 2--59
## [61] 2--61 3--5 3--6 3--8 3--10 3--11 3--13 3--14 3--16 3--17 3--19 3--20
## [73] 3--21 3--24 3--25 3--26 3--27 3--28 3--29 3--30 3--32 3--33 3--34 3--35
## [85] 3--36 3--37 3--38 3--39 3--41 3--42 3--44 3--45 3--46 3--48 3--49 3--51
## [97] 3--53 3--54 3--55 3--56 3--57 3--59 3--60 3--62 4--5 4--8 4--9 4--10
## [109] 4--12 4--14 4--15 4--16 4--17 4--21 4--25 4--29 4--30 4--34 4--35 4--39
## + ... omitted several edges
par(mfrow=c(1,2), mar = c(2, 2, 2, 2))
kc <- fastgreedy.community(mmd_graph)
sizes(kc)
## Community sizes
## 1 2
## 32 30
plot(kc, mmd_graph, vertex.labels=NA, mar = c(8, 8, 8, 8))
dendPlot(kc, mode="phylo")
# Applying the Fast Greedy algorithm for community detection
kc <- fastgreedy.community(mmd_graph)
# Getting the community membership for each vertex
community_membership <- as.vector(membership(kc))
# Create a data frame with the vertices and their corresponding community membership
vertex_community_df <- data.frame(vertex = V(mmd_graph)$name, community = community_membership)
# Print the vertex-community mapping
print(vertex_community_df)
## vertex community
## 1 1 2
## 2 2 2
## 3 3 1
## 4 4 2
## 5 5 2
## 6 6 2
## 7 7 2
## 8 8 2
## 9 9 2
## 10 10 2
## 11 11 1
## 12 12 2
## 13 13 2
## 14 14 2
## 15 15 1
## 16 16 1
## 17 17 2
## 18 18 2
## 19 19 2
## 20 20 1
## 21 21 2
## 22 22 1
## 23 23 2
## 24 24 2
## 25 25 1
## 26 26 1
## 27 27 1
## 28 28 1
## 29 29 1
## 30 30 1
## 31 31 1
## 32 32 1
## 33 33 1
## 34 34 1
## 35 35 1
## 36 36 1
## 37 37 2
## 38 38 2
## 39 39 2
## 40 40 1
## 41 41 1
## 42 42 1
## 43 43 1
## 44 44 1
## 45 45 1
## 46 46 1
## 47 47 2
## 48 48 1
## 49 49 1
## 50 50 2
## 51 51 2
## 52 52 2
## 53 53 2
## 54 54 1
## 55 55 2
## 56 56 1
## 57 57 1
## 58 58 2
## 59 59 1
## 60 60 1
## 61 61 2
## 62 62 2
# Change the values 2 and 1 to 0 and 1, respectively
vertex_community_df$community <- ifelse(vertex_community_df$community == 2, 1, 0)
# Print the modified vertex-community mapping
print(vertex_community_df)
## vertex community
## 1 1 1
## 2 2 1
## 3 3 0
## 4 4 1
## 5 5 1
## 6 6 1
## 7 7 1
## 8 8 1
## 9 9 1
## 10 10 1
## 11 11 0
## 12 12 1
## 13 13 1
## 14 14 1
## 15 15 0
## 16 16 0
## 17 17 1
## 18 18 1
## 19 19 1
## 20 20 0
## 21 21 1
## 22 22 0
## 23 23 1
## 24 24 1
## 25 25 0
## 26 26 0
## 27 27 0
## 28 28 0
## 29 29 0
## 30 30 0
## 31 31 0
## 32 32 0
## 33 33 0
## 34 34 0
## 35 35 0
## 36 36 0
## 37 37 1
## 38 38 1
## 39 39 1
## 40 40 0
## 41 41 0
## 42 42 0
## 43 43 0
## 44 44 0
## 45 45 0
## 46 46 0
## 47 47 1
## 48 48 0
## 49 49 0
## 50 50 1
## 51 51 1
## 52 52 1
## 53 53 1
## 54 54 0
## 55 55 1
## 56 56 0
## 57 57 0
## 58 58 1
## 59 59 0
## 60 60 0
## 61 61 1
## 62 62 1
Setting a seed:
set.seed(5114382)
Here, we started with a plot reporting the full network of our data for having a general idea: it’s possible to see that there are many edges and so many connections between the nodes, below we’ll analyze deeply these connections in terms of betweenness, closeness and other parameters of the nodes.
plot(mmd_graph, layout= layout_with_fr(mmd_graph),
vertex.color = "lightpink",
edge.color = "lightblue",
main = "Network Plot",
edge.curved = 0.1
)
A way to combine these ‘multiedges’ into weighted edges. That is, each pair of individuals will be connected by one edge, but the number of interactions will be represented by the edge widths.
E(mmd_graph)$weight=1 #make all edge weights = 1
ig2=simplify(mmd_graph, remove.multiple=T, edge.attr.comb=list(weight=sum))
plot(ig2, edge.width=E(ig2)$weight,
vertex.color = "lightblue")
Then we proceeded with many plots applying a different layout each time, with the aim of visualize the full network from different point of view.
par(mfrow=c(1, 2))
# Define colors for nodes
node_colors <- c("red", "blue", "green", "yellow", "orange")
vertex_border_colors = "blue"
vertex_size <- 15
# Plot the graph with Fruchterman-Reingold layout
plot(mmd_graph, layout = layout_with_fr(mmd_graph),
vertex.color = "lightblue",
vertex.frame.color = node_colors[V(mmd_graph)$color], # fill nodes with colors
vertex.size = vertex_size,
vertex.label.color = "blue",
vertex.label.dist = 0.5,
vertex.label.cex = 0.8,
vertex.border.color = vertex_border_colors,
edge.width = E(mmd_graph)$width,
edge.color = "purple",
edge.curved = 2.0, # Adjust the edge distance (0.2 is just an example)
main = "Network Plot"
)
# Create the Fruchterman-Reingold layout
layout <- layout.fruchterman.reingold(mmd_graph)
# Set the desired colors for the nodes
node_colors <- c("purple", "lightblue", "blue", "yellow")
# Plot the graph with the Fruchterman-Reingold layout and custom node colors
plot(mmd_graph, layout = layout, vertex.label = NA, vertex.color = node_colors,
edge_colors = "lightpink")
library(igraph)
library(plotly)
## Warning: il pacchetto 'plotly' è stato creato con R versione 4.2.3
##
## Caricamento pacchetto: 'plotly'
## Il seguente oggetto è mascherato da 'package:igraph':
##
## groups
## Il seguente oggetto è mascherato da 'package:ggplot2':
##
## last_plot
## Il seguente oggetto è mascherato da 'package:stats':
##
## filter
## Il seguente oggetto è mascherato da 'package:graphics':
##
## layout
# Create the Fruchterman-Reingold layout
layout <- layout.fruchterman.reingold(mmd_graph)
# Create an interactive 3D plot
plot_ly(
type = "scatter3d",
x = layout[, 1],
y = layout[, 2],
z = layout[, 2],
mode = "markers",
marker = list(
size = 5,
color = "purple",
line = list(color = "lightblue", width = 0.5)
),
hoverinfo = "text",
text = V(mmd_graph)$name
) %>%
layout(
scene = list(
xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Z")
)
)
vertex_size <- 10
plot(mmd_graph, layout = layout_with_fr(mmd_graph),
vertex.color = "lightblue",
vertex.frame.color = node_colors[V(mmd_graph)$color], # Fill nodes with colors
vertex.size = vertex_size,
vertex.label.color = "blue",
vertex.label.dist = 0.8,
vertex.label.cex = 0.8,
vertex.border.color = vertex_border_colors,
edge.width = E(mmd_graph)$width,
edge.color = "purple",
main = "Network Plot"
)
# Create the circular layout
layout <- layout_in_circle(mmd_graph)
# Plot the graph with the circular layout
node_colors <- c("purple", "lightblue", "blue", "yellow")
inside_color <- "lightblue"
plot(mmd_graph, layout = layout, vertex.label = NA, vertex.color = node_colors,
vertex.frame.color = inside_color)
library(plotly)
# Create the circular layout
layout <- layout_in_circle(mmd_graph)
# Create an interactive 3D plot
plot_ly(
type = "scatter3d",
x = layout[, 1],
y = layout[, 2],
z = layout[, 2],
mode = "markers",
marker = list(
size = 5,
color = "purple",
line = list(color = "black", width = 0.5)
),
hoverinfo = "text",
text = V(mmd_graph)$name
) %>%
layout(
scene = list(
xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Z")
)
)
# Create the Kamada-Kawai layout
layout <- layout.kamada.kawai(mmd_graph)
# Plot the graph with the Kamada-Kawai layout
node_colors <- c("purple", "lightblue", "blue", "yellow")
plot(mmd_graph, layout = layout, vertex.label = NA, vertex.color = node_colors)
mmd_degree <- igraph::degree(mmd_graph)
degree_distribution(mmd_graph)
## [1] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [7] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [13] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.01612903
## [19] 0.00000000 0.00000000 0.01612903 0.00000000 0.00000000 0.00000000
## [25] 0.01612903 0.00000000 0.00000000 0.01612903 0.01612903 0.03225806
## [31] 0.01612903 0.04838710 0.04838710 0.01612903 0.09677419 0.08064516
## [37] 0.04838710 0.06451613 0.06451613 0.03225806 0.03225806 0.01612903
## [43] 0.01612903 0.04838710 0.06451613 0.04838710 0.04838710 0.03225806
## [49] 0.00000000 0.01612903 0.00000000 0.01612903 0.00000000 0.01612903
## [55] 0.00000000 0.01612903
sort(mmd_degree, decreasing = TRUE)[1:5]
## 48 30 59 12 16
## 55 53 51 49 47
It seems that the 5 nodes returned as outputs are the most important in terms of information flow, influence and in some cases, like in this one, it can represents an index of the spreadness of a desease.
par(mfrow=c(1, 2))
plot(mmd_degree, V(mmd_graph), pch = 15)
points(mmd_degree, V(mmd_graph), col = "purple", pch = 15)
plot(mmd_degree, V(mmd_graph), type = "o", col = "blue", pch = 16)
points(mmd_degree, V(mmd_graph), col = "purple", pch = 18)
igraph.options(mmd_graph, default = NULL )
netstat_dist <- function (s, ...) {
h <- hist(s, plot = FALSE, ...)
bin_breaks <- as.vector(stats::filter(h$breaks, c(1, 1) / 2))
bin_breaks <- bin_breaks[-length(h$breaks)]
list(s = bin_breaks[h$counts > 0], fs = h$counts[h$counts > 0])
}
par(mfrow = c(1, 3))
hist(degree(mmd_graph), col = "lightblue")
hist(graph.strength(mmd_graph), col = "lightpink")
The distributions are perfectly the same, indeed they follow a gaussian distribution. This fact suggests that the degree are homgenously distributed meaning that the nodes have the same number of connections. moreover, the strength of a node’s connection is not influenced by its degree.
ks_test <- ks.test(degree(mmd_graph), graph.strength(mmd_graph))
ks_test
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: degree(mmd_graph) and graph.strength(mmd_graph)
## D = 0, p-value = 1
## alternative hypothesis: two-sided
We also perfomr the KS test i order to have a more impact proof in order to show that we are dealing with two equal distributions. Indeed, the p-value is above 0.5.
par(mfrow = c(2, 1))
# Compute the degree distribution
degree_dist <- degree(mmd_graph)
# Set the colors for the bars
bar_colors <- ifelse(degree_dist >= 10, "lightblue", "gray")
# Plot the degree distribution histogram with colors
barplot(degree_dist, col = bar_colors, main = "Degree Distribution", xlab = "Degree", ylab = "Frequency")
# Compute the node strength distribution
node_strength <- graph.strength(mmd_graph)
# Set the colors for the bars
bar_colors <- ifelse(node_strength >= 50, "lightpink", "lightpink")
# Plot the node strength distribution histogram with colors
barplot(node_strength, col = bar_colors, main = "Node Strength Distribution", xlab = "Node", ylab = "Strength")
# Generate the degree sequence
degrees <- degree(mmd_graph)
# Set up the color vector
colors <- ifelse(degrees > 0, "blue", "red") # Color nodes with degree > 0 as blue, and degree 0 as red
# Plot the degree distribution with colors
plot(degrees, pch = 16, col = colors, main = "Degree Distribution", xlab = "Node", ylab = "Degree")
mmd_close <- igraph::closeness(mmd_graph, mode="total")
sort(mmd_close, decreasing=TRUE)[1:5]
## 48 30 59 12 16
## 0.01492537 0.01449275 0.01408451 0.01369863 0.01333333
The first five nodes have a quick interaction with all the others since they have a shorter communication path.
par(mfrow=c(1,2))
plot(mmd_close, V(mmd_graph), col= "blue", pch = 18)
points(mmd_degree, V(mmd_graph), col = "purple", pch = 18)
plot(mmd_close, V(mmd_graph), type = "o", col = "blue", pch = 16)
points(mmd_degree, V(mmd_graph), col = "purple", pch = 18)
mmd_betweeness <- igraph::betweenness(mmd_graph)
sorted_betweenness <- sort(mmd_betweeness, decreasing = TRUE)
sorted_betweenness
## 48 30 59 12 55 13 40 14
## 30.319064 25.465322 23.022910 21.799489 20.557544 19.549064 19.232791 19.021568
## 62 16 56 53 39 3 25 61
## 18.060715 17.130029 16.874347 16.410616 16.007716 15.824153 15.480694 15.368677
## 49 45 41 24 35 10 36 11
## 15.348905 14.854001 14.668664 14.362922 13.568335 13.166998 11.848321 11.605552
## 51 8 31 42 19 33 37 23
## 11.553339 11.532951 11.471701 11.317816 11.243507 10.931635 10.495686 10.352761
## 5 44 50 1 60 20 58 57
## 10.201680 10.075797 9.638237 9.451681 9.348316 9.188040 9.072970 8.802223
## 29 47 54 38 9 52 26 17
## 8.778183 8.701615 8.437928 8.382912 8.041797 7.917197 7.393796 7.390721
## 43 27 21 34 18 22 28 32
## 7.375403 7.295375 7.149872 6.912594 6.835462 6.808295 5.949095 5.332225
## 15 4 46 2 7 6
## 5.323145 4.473049 4.442393 3.078531 2.004435 1.749241
names(sorted_betweenness)[1:5]
## [1] "48" "30" "59" "12" "55"
These 5 nodes have the greatest number of shortest paths between the other nodes; they are the most important connectors. Plotting the betweenness centrality we can state that all the pints inside the graph are blue, it means that the values are all grater than 0 and they indicates their importance in connecting the different part of the network.
s <- graph.strength(mmd_graph)
b <- betweenness(mmd_graph)
# Set up the color vector
colors <- ifelse(b > 0, "blue", "red") # Color nodes with betweenness > 0 as blue, and betweenness 0 as red
# Plot graph strength against betweenness centrality with colors
plot(s, b, pch = 16, col = colors, main = "Graph Strength vs Betweenness Centrality", xlab = "Graph Strength", ylab = "Betweenness")
The plot shows a remarkable relationship between the graph stregth and the betweennss, more precisely they have a positive correlation. The nodes that have both high betweenness and streght are important hubs since they play a crucial role in facilitating communication and information flows.
transitivity(mmd_graph, type="global")
## [1] 0.659942
transitivity(mmd_graph, type="local")[1:5]
## [1] 0.6884780 0.7318841 0.6553911 0.7389163 0.6327986
The obtained results suggest that exists a high level of cohesion among subgroups and communities.
mmd_pr <- igraph::page.rank(mmd_graph)$vector
sort(mmd_pr, decreasing=TRUE)[1:4]
## 48 30 59 12
## 0.02257253 0.02178759 0.02104352 0.02033998
It highlights the nodes with the highest centrality and influence, in this case the most dominant ones based on their connectivity.
igraph::diameter(mmd_graph)
## [1] 2
The greatest number of steps required to reach any node inside the network.
mean_d <- mean_distance(mmd_graph)
mean_d
## [1] 1.382866
D <- distances(mmd_graph)
D
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## 1 0 1 1 2 2 2 2 1 2 1 1 1 1 1 2 1 2 1 2 1 1 1 1 1 2 1 2 2
## 2 1 0 2 2 2 2 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
## 3 1 2 0 2 1 1 2 1 2 1 1 2 1 1 2 1 1 2 1 1 1 2 2 1 1 1 1 1
## 4 2 2 2 0 1 2 2 1 1 1 2 1 2 1 1 1 1 2 2 2 1 2 2 2 1 2 2 2
## 5 2 2 1 1 0 2 1 1 1 1 2 2 1 2 1 1 2 1 2 2 1 2 2 2 1 1 2 2
## 6 2 2 1 2 2 0 2 2 2 1 2 1 1 1 2 2 2 1 2 2 2 2 2 1 2 2 2 2
## 7 2 2 2 2 1 2 0 1 2 1 2 2 1 2 2 2 1 2 1 2 2 2 1 2 1 2 2 2
## 8 1 2 1 1 1 2 1 0 1 1 2 2 2 1 2 1 1 1 1 1 1 2 2 2 1 2 2 1
## 9 2 2 2 1 1 2 2 1 0 1 1 1 1 2 1 1 1 2 1 1 1 2 2 1 1 2 2 2
## 10 1 1 1 1 1 1 1 1 1 0 2 1 1 1 2 2 1 2 1 2 1 1 2 1 1 1 2 2
## 11 1 2 1 2 2 2 2 2 1 2 0 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 2
## 12 1 1 2 1 2 1 2 2 1 1 1 0 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## 13 1 2 1 2 1 1 1 2 1 1 1 2 0 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2
## 14 1 2 1 1 2 1 2 1 2 1 1 1 1 0 2 1 1 2 1 1 1 1 1 1 1 1 1 1
## 15 2 2 2 1 1 2 2 2 1 2 2 1 2 2 0 1 1 2 1 1 2 1 1 2 1 2 1 2
## 16 1 2 1 1 1 2 2 1 1 2 1 1 1 1 1 0 1 2 1 1 1 1 1 1 1 2 1 1
## 17 2 2 1 1 2 2 1 1 1 1 2 2 2 1 1 1 0 2 2 2 1 2 1 2 1 1 2 1
## 18 1 2 2 2 1 1 2 1 2 2 1 1 1 2 2 2 2 0 2 2 2 2 2 1 2 2 1 1
## 19 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 2 2 0 2 1 1 1 1 1 2 2 1
## 20 1 2 1 2 2 2 2 1 1 2 2 1 1 1 1 1 2 2 2 0 1 1 1 1 1 2 2 1
## 21 1 2 1 1 1 2 2 1 1 1 1 1 1 1 2 1 1 2 1 1 0 1 1 1 1 2 2 2
## 22 1 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 2 1 1 1 0 1 1 1 1 2 2
## 23 1 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 1 2 1 1 1 1 0 1 1 1 1 1
## 24 1 1 1 2 2 1 2 2 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 0 2 2 1 1
## 25 2 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 0 1 2 1
## 26 1 2 1 2 1 2 2 2 2 1 1 1 1 1 2 2 1 2 2 2 2 1 1 2 1 0 2 2
## 27 2 2 1 2 2 2 2 2 2 2 1 1 1 1 1 1 2 1 2 2 2 2 1 1 2 2 0 1
## 28 2 2 1 2 2 2 2 1 2 2 2 1 2 1 2 1 1 1 1 1 2 2 1 1 1 2 1 0
## 29 2 1 1 1 2 2 2 2 2 2 1 2 2 2 1 1 2 1 1 1 1 1 2 2 1 1 1 2
## 30 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1
## 31 1 2 2 2 2 2 2 1 2 2 2 1 1 1 2 1 2 1 1 2 2 1 1 1 1 1 1 1
## 32 2 2 1 2 1 2 2 1 2 2 1 2 1 1 2 1 2 1 2 1 2 2 2 2 1 2 1 1
## 33 1 2 1 2 2 1 1 1 2 2 1 1 1 2 2 1 1 2 2 2 2 2 2 1 1 1 1 1
## 34 2 2 1 1 2 2 2 2 1 2 1 2 2 1 2 2 2 2 1 2 2 1 2 1 1 1 1 2
## 35 1 1 1 1 2 2 2 1 2 2 1 1 1 1 2 1 2 1 2 1 2 2 1 1 1 1 1 1
## 36 1 1 1 2 2 2 2 2 2 1 1 1 2 1 1 1 2 1 2 1 1 2 1 1 2 2 1 1
## 37 1 1 1 2 1 1 2 2 2 1 2 1 2 1 2 2 1 2 1 2 2 1 1 1 2 1 2 2
## 38 1 1 1 2 1 2 1 1 2 2 2 1 1 1 2 1 2 1 2 2 2 2 2 1 2 2 2 1
## 39 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 2 2 2
## 40 1 1 2 2 2 2 2 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 1 1 1 1
## 41 2 1 1 2 1 2 2 2 2 1 1 1 1 2 2 1 2 1 1 1 1 1 1 1 1 1 1 1
## 42 1 2 1 2 1 2 2 2 1 2 2 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 2
## 43 1 2 2 2 1 2 2 2 2 2 1 1 1 1 2 1 2 2 1 1 2 1 2 1 1 2 1 2
## 44 2 1 1 1 2 2 2 2 1 2 2 1 1 2 1 1 2 1 2 1 2 1 2 2 2 1 1 1
## 45 1 1 1 1 1 2 2 2 1 1 1 1 1 1 2 1 2 1 2 1 1 1 2 2 1 1 1 2
## 46 2 2 1 2 2 2 2 2 1 2 2 1 1 1 1 2 2 2 2 1 2 2 2 1 1 2 2 2
## 47 1 1 2 2 1 2 2 1 1 2 1 1 1 1 2 1 1 1 2 2 2 2 1 1 1 2 1 1
## 48 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## 49 1 2 1 1 1 2 1 2 1 1 1 2 1 1 1 1 1 2 1 1 2 2 2 2 1 1 1 2
## 50 1 1 2 2 1 1 2 1 1 2 1 1 1 1 2 1 2 1 1 1 2 1 1 2 2 1 1 2
## 51 1 1 1 2 1 2 2 1 2 1 1 1 1 1 2 2 1 1 1 1 2 1 1 1 2 2 1 2
## 52 1 2 2 2 2 1 2 1 2 1 1 1 1 1 2 2 1 2 1 2 1 2 1 1 2 1 1 1
## 53 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 2 1 2 1 1 2 1 1 1 2 2
## 54 2 2 1 1 2 2 2 2 1 1 2 1 1 1 1 1 2 2 1 1 1 1 2 2 1 1 2 1
## 55 1 1 1 1 1 2 1 1 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 2 1 1 1
## 56 1 1 1 1 1 2 1 1 1 1 2 1 2 2 1 1 1 2 1 1 2 1 2 1 1 1 1 1
## 57 1 2 1 2 2 2 2 1 2 1 1 1 1 2 2 1 2 1 2 1 2 2 1 1 2 2 1 1
## 58 1 1 2 1 1 2 1 2 1 1 1 2 1 1 2 2 1 2 1 2 1 2 1 2 2 2 2 2
## 59 1 1 1 1 1 2 1 1 2 2 1 1 1 1 2 2 2 1 1 1 2 1 1 1 1 1 2 1
## 60 2 2 1 1 1 2 1 1 1 1 1 2 2 1 1 1 1 2 2 2 2 2 2 2 1 1 2 2
## 61 2 1 2 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 1 2 1 1 2 2 2 2 2
## 62 2 2 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 2 1 1 2 2
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## 1 2 1 1 2 1 2 1 1 1 1 2 1 2 1 1 2 1 2 1 1 1 1 1 1 2
## 2 1 1 2 2 2 2 1 1 1 1 2 1 1 2 2 1 1 2 1 1 2 1 1 2 1
## 3 1 1 2 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 2 1 1 2 1 2 1
## 4 1 1 2 2 2 1 1 2 2 2 1 2 2 2 2 1 1 2 2 1 1 2 2 2 1
## 5 2 1 2 1 2 2 2 2 1 1 1 2 1 1 1 2 1 2 1 2 1 1 1 2 1
## 6 2 1 2 2 1 2 2 2 1 2 1 2 2 2 2 2 2 2 2 1 2 1 2 1 1
## 7 2 2 2 2 1 2 2 2 2 1 1 2 2 2 2 2 2 2 2 1 1 2 2 2 2
## 8 2 1 1 1 1 2 1 2 2 1 1 2 2 2 2 2 2 2 1 1 2 1 1 1 1
## 9 2 2 2 2 2 1 2 2 2 2 1 1 2 1 2 1 1 1 1 1 1 1 2 2 1
## 10 2 1 2 2 2 2 2 1 1 2 1 1 1 2 2 2 1 2 2 2 1 2 1 1 1
## 11 1 1 2 1 1 1 1 1 2 2 1 1 1 2 1 2 1 2 1 1 1 1 1 1 1
## 12 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
## 13 2 1 1 1 1 2 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 14 2 1 1 1 2 1 1 1 1 1 1 2 2 1 1 2 1 1 1 1 1 1 1 1 1
## 15 1 1 2 2 2 2 2 1 2 2 1 1 2 1 2 1 2 1 2 1 1 2 2 2 1
## 16 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 2 1 1 1 1 2 2 1
## 17 2 2 2 2 1 2 2 2 1 2 1 1 2 1 2 2 2 2 1 2 1 2 1 1 2
## 18 1 1 1 1 2 2 1 1 2 1 1 1 1 1 2 1 1 2 1 1 2 1 1 2 2
## 19 1 2 1 2 2 1 2 2 1 2 1 1 1 2 1 2 2 2 2 1 1 1 1 1 1
## 20 1 1 2 1 2 2 1 1 2 2 2 1 1 1 1 1 1 1 2 1 1 1 1 2 2
## 21 1 1 2 2 2 2 2 1 2 2 1 1 1 1 2 2 1 2 2 1 2 2 2 1 1
## 22 1 1 1 2 2 1 2 2 1 2 2 1 1 2 1 1 1 2 2 1 2 1 1 2 1
## 23 2 1 1 2 2 2 1 1 1 2 2 2 1 1 2 2 2 2 1 1 2 1 1 1 2
## 24 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 2 1 1 1
## 25 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 2 1 1 1 1 1 2 2 2 1
## 26 1 1 1 2 1 1 1 2 1 2 2 1 1 1 2 1 1 2 2 1 1 1 2 1 1
## 27 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 2 1 1 1 1 1 1 2
## 28 2 1 1 1 1 2 1 1 2 1 2 1 1 2 2 1 2 2 1 1 2 2 2 1 2
## 29 0 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 2 2 1 1 2 1 2 2
## 30 1 0 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1
## 31 1 1 0 1 1 1 1 1 1 2 2 1 1 1 1 1 2 1 1 1 2 2 1 2 1
## 32 1 1 1 0 1 1 1 1 2 2 2 1 2 2 1 2 1 1 2 1 2 2 2 2 1
## 33 1 2 1 1 0 1 1 2 1 1 2 1 1 1 1 1 1 2 2 1 1 2 2 2 2
## 34 1 1 1 1 1 0 1 2 2 1 1 1 1 1 1 1 1 2 2 1 1 2 2 2 2
## 35 1 1 1 1 1 1 0 1 2 1 1 1 1 1 1 1 1 2 2 1 1 2 1 2 1
## 36 1 1 1 1 2 2 1 0 1 1 2 2 1 2 1 1 1 1 2 1 1 1 1 1 1
## 37 1 1 1 2 1 2 2 1 0 1 2 1 2 1 2 1 1 2 1 1 1 2 1 1 1
## 38 1 1 2 2 1 1 1 1 1 0 1 1 2 2 2 2 2 2 1 1 2 1 1 1 2
## 39 2 1 2 2 2 1 1 2 2 1 0 1 1 1 1 1 1 2 2 1 1 1 1 1 1
## 40 1 1 1 1 1 1 1 2 1 1 1 0 1 2 1 1 1 1 1 1 1 1 2 2 1
## 41 1 1 1 2 1 1 1 1 2 2 1 1 0 1 1 1 1 1 2 1 1 1 1 2 1
## 42 2 1 1 2 1 1 1 2 1 2 1 2 1 0 1 1 1 2 2 1 1 1 1 2 1
## 43 1 1 1 1 1 1 1 1 2 2 1 1 1 1 0 2 1 2 1 1 2 1 2 1 2
## 44 1 1 1 2 1 1 1 1 1 2 1 1 1 1 2 0 1 1 1 1 1 2 1 2 1
## 45 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 0 1 2 1 2 2 1 2 1
## 46 2 1 1 1 2 2 2 1 2 2 2 1 1 2 2 1 1 0 2 2 1 1 2 1 1
## 47 2 2 1 2 2 2 2 2 1 1 2 1 2 2 1 1 2 2 0 1 1 1 1 1 1
## 48 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 0 1 1 1 1 1
## 49 1 1 2 2 1 1 1 1 1 2 1 1 1 1 2 1 2 1 1 1 0 1 1 2 1
## 50 2 2 2 2 2 2 2 1 2 1 1 1 1 1 1 2 2 1 1 1 1 0 1 2 2
## 51 1 1 1 2 2 2 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 0 1 1
## 52 2 1 2 2 2 2 2 1 1 1 1 2 2 2 1 2 2 1 1 1 2 2 1 0 1
## 53 2 1 1 1 2 2 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 2 1 1 0
## 54 1 1 1 1 2 1 1 1 2 2 2 1 1 1 1 2 1 2 2 1 2 1 2 2 1
## 55 2 1 1 2 1 2 1 2 2 1 1 1 1 2 1 2 2 1 2 2 1 1 1 1 1
## 56 1 1 2 2 1 1 1 1 2 1 1 2 1 2 2 1 1 1 2 1 1 2 2 2 1
## 57 2 1 2 1 1 1 1 1 2 1 2 1 1 1 1 1 1 2 1 1 2 2 1 2 2
## 58 2 1 1 1 2 2 1 2 1 2 1 2 2 2 2 2 1 2 2 1 1 2 1 1 2
## 59 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1
## 60 2 1 1 2 1 1 1 2 1 2 1 1 1 1 2 2 1 1 2 1 1 2 2 2 1
## 61 2 2 1 2 2 1 2 2 1 1 1 1 1 1 1 2 1 1 1 2 1 2 2 2 1
## 62 2 1 1 1 2 1 1 1 1 2 1 1 1 1 2 1 2 1 1 1 1 2 2 1 1
## 54 55 56 57 58 59 60 61 62
## 1 2 1 1 1 1 1 2 2 2
## 2 2 1 1 2 1 1 2 1 2
## 3 1 1 1 1 2 1 1 2 1
## 4 1 1 1 2 1 1 1 1 1
## 5 2 1 1 2 1 1 1 1 1
## 6 2 2 2 2 2 2 2 1 1
## 7 2 1 1 2 1 1 1 1 1
## 8 2 1 1 1 2 1 1 1 1
## 9 1 1 1 2 1 2 1 1 1
## 10 1 1 1 1 1 2 1 1 1
## 11 2 2 2 1 1 1 1 1 1
## 12 1 1 1 1 2 1 2 1 1
## 13 1 1 2 1 1 1 2 2 2
## 14 1 2 2 2 1 1 1 2 1
## 15 1 1 1 2 2 2 1 1 2
## 16 1 2 1 1 2 2 1 1 1
## 17 2 1 1 2 1 2 1 1 1
## 18 2 1 2 1 2 1 2 2 1
## 19 1 1 1 2 1 1 2 1 1
## 20 1 1 1 1 2 1 2 1 1
## 21 1 1 2 2 1 2 2 2 1
## 22 1 1 1 2 2 1 2 1 1
## 23 2 1 2 1 1 1 2 1 1
## 24 2 1 1 1 2 1 2 2 2
## 25 1 2 1 2 2 1 1 2 1
## 26 1 1 1 2 2 1 1 2 1
## 27 2 1 1 1 2 2 2 2 2
## 28 1 1 1 1 2 1 2 2 2
## 29 1 2 1 2 2 1 2 2 2
## 30 1 1 1 1 1 1 1 2 1
## 31 1 1 2 2 1 1 1 1 1
## 32 1 2 2 1 1 1 2 2 1
## 33 2 1 1 1 2 1 1 2 2
## 34 1 2 1 1 2 1 1 1 1
## 35 1 1 1 1 1 1 1 2 1
## 36 1 2 1 1 2 1 2 2 1
## 37 2 2 2 2 1 1 1 1 1
## 38 2 1 1 1 2 1 2 1 2
## 39 2 1 1 2 1 1 1 1 1
## 40 1 1 2 1 2 1 1 1 1
## 41 1 1 1 1 2 1 1 1 1
## 42 1 2 2 1 2 1 1 1 1
## 43 1 1 2 1 2 1 2 1 2
## 44 2 2 1 1 2 1 2 2 1
## 45 1 2 1 1 1 1 1 1 2
## 46 2 1 1 2 2 1 1 1 1
## 47 2 2 2 1 2 1 2 1 1
## 48 1 2 1 1 1 1 1 2 1
## 49 2 1 1 2 1 1 1 1 1
## 50 1 1 2 2 2 2 2 2 2
## 51 2 1 2 1 1 2 2 2 2
## 52 2 1 2 2 1 1 2 2 1
## 53 1 1 1 2 2 1 1 1 1
## 54 0 2 2 1 1 1 1 2 1
## 55 2 0 1 1 1 1 1 1 1
## 56 2 1 0 1 1 1 1 1 1
## 57 1 1 1 0 1 1 1 1 2
## 58 1 1 1 1 0 1 1 1 2
## 59 1 1 1 1 1 0 1 1 1
## 60 1 1 1 1 1 1 0 1 1
## 61 2 1 1 1 1 1 1 0 2
## 62 1 1 1 2 2 1 1 2 0
Minimum numbers of edges required to travel from one node to another.
igraph::max_cliques(mmd_graph)%>% map_dbl(length)%>%table()
## .
## 5 6 7 8 9 10 11 12 13
## 15 119 482 848 1273 1480 1033 328 50
It shows the cliques table in which there is the number of times each length occurs.
Let’s check which is the most important node standing for the most important, so dominant individual:
ia <- order(b, decreasing = T)[1] # we order from maximum to minimum in order to get the most important
V(mmd_graph)$name[ia]
## [1] "48"
g_raph <- subgraph.edges(mmd_graph, E(mmd_graph)[.inc(ia)])
plot(g_raph, vertex.label = NA, vertex.color = c("purple", "lightblue", "lightgreen", "yellow"))
library(ergm)
am <- get.adjacency(mmd_graph, sparse = FALSE)
g <- as.network(am, directed= FALSE)
ergm(g~edges)%>% summary()
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Evaluating log-likelihood at the estimate.
## Call:
## ergm(formula = g ~ edges)
##
## Maximum Likelihood Results:
##
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges 0.47740 0.04731 0 10.09 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 2621 on 1891 degrees of freedom
## Residual Deviance: 2517 on 1890 degrees of freedom
##
## AIC: 2519 BIC: 2524 (Smaller is better. MC Std. Err. = 0)
Since the degrees are homogenous we have fitted the ERG model. Based on the output there is a positive effect of the edges on the likelihood. The model fits the data better than the null one.
degree_dist <- function (mmd_graph) {
fd <- table(degree(mmd_graph))
d <- as.numeric(names(fd)) + 1 # degree + 1
list(d = d, fd = fd)
}
dd <- degree_dist(mmd_graph)
with(dd, plot(log(d), log(fd)))
We have performed a scatter plot of the degrees of the logarithmic scale in order to better understand the distribution’s shape. It appears that the degree follows a poisson distribution, for that reason below we fit two models, the naive linear model and the generalized linear model specifying the poisson family.
dd <- degree_dist(mmd_graph)
m0 <- lm(log(fd)~log(d), data = dd) # naive linear regression
m1 <- glm(fd~log(d), family = poisson, data = dd) # poisson regression
summary(m0)
##
## Call:
## lm(formula = log(fd) ~ log(d), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75984 -0.54598 0.02098 0.50550 1.20880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7550 1.5769 -0.479 0.636
## log(d) 0.3763 0.4353 0.865 0.395
##
## Residual standard error: 0.6265 on 26 degrees of freedom
## Multiple R-squared: 0.02794, Adjusted R-squared: -0.009444
## F-statistic: 0.7474 on 1 and 26 DF, p-value: 0.3952
summary(m1)
##
## Call:
## glm(formula = fd ~ log(d), family = poisson, data = dd)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0826 -0.8490 -0.2228 0.5203 2.1329
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3133 1.7568 -0.178 0.858
## log(d) 0.3058 0.4821 0.634 0.526
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 22.936 on 27 degrees of freedom
## Residual deviance: 22.524 on 26 degrees of freedom
## AIC: 97.866
##
## Number of Fisher Scoring iterations: 5
# Perform ANOVA to compare the models
anova_result <- anova(m0, m1)
print(anova_result)
## Analysis of Variance Table
##
## Response: log(fd)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(d) 1 0.2934 0.29340 0.7474 0.3952
## Residuals 26 10.2065 0.39256